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We present a hybrid quantum-classical algorithm to simulate thermal states of classical Hamil- 
tonians on a quantum computer. Our scheme employs a sequence of locally controlled rotations, 
building up the desired state by adding qubits one at a time. We identify a class of classical models 
for which our method is efficient and avoids potential exponential overheads encountered by Grover- 
like or quantum Metropolis schemes. Our algorithm also gives an exponential advantage for 2D 
Ising models with magnetic field on a square lattice, compared with the previously known Zalka's 
algorithm. 
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Simulation of a finite-temperature physical system 
with a controllable quantum device is one of the most 
important goals of quantum simulation [Jj |2j. Classi- 
cal Markov-Chain Monte Carlo (MCMC) algorithms are 
powerful tools for sampling Gibbs distributions. They 
are efficient provided that the gap A of the transition 
matrix is non- vanishing; the running time typically scales 
as r ~ O (1/^)- A quantum generalization [3 of MCMC 
has recently been explored by the quantum information 
community due to the connection to quantum walks [4]. 
Richter [5 developed a method for sampling from the 
Gibbs distribution for periodic lattices. Somma et at. 
[6^ combined quantum walk and quantum Zeno effect to 
achieve quantum speedup. Wocjan and Abeyesinghe [7] 
improved it by using fixed point quantum search. Gen- 
erally, these quantum algorithms allow the running time 
to scale as r ~ 0(1/V^), a quadratic speedup compared 
with the classical counterparts. However, for many prob- 
lems of practical interest, such as optimization problems 
and spin glasses, the gap A may become exponentially 
small when the system size increases, making it unprac- 
tical to use MCMC algorithms for solving them (see Fig. 
fl]). Therefore, gap-independent methods are more desir- 
able for solving these problems. 

A class of gap-independent methods is called belief 
propagation [8], which generalizes the transfer matrix 
methods in statistical physics. For problems involving 
a regular geometry, it can be very efficient, e.g. one- 
dimensional spin chains. This property will be exploited 
in this letter, where a different way for obtaining samples 
from the thermal state is discussed. This approach is a 
generalization of the state preparation method by Lidar 
and Biham [9 , and Zalka [10]. We show that in some 
cases, the structure of the system under investigation al- 
lows for large speedups over the general methods. This 
is because the cost of our method is independent of the 
temperature and the gap size. 
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A* Markov Chain Gap A 

FIG. 1: (color online) The running time r ^ O (1/A) of 
Markov chain methods is limited by the gap A of the Markov 
matrix. A quantum quadratic speedup r ^ 0(1/a/A) (black 
solid line) relative to classical Markov chains (black dashed 
line) can be achieved by a quantum computer. Below some 
critical gap size A < A*, Markov chain methods become inef- 
ficient (shaded region), and classical belief propagation meth- 
ods (red dashed line), including transfer matrix methods, be- 
come more efficient. Combined with amplitude amplification, 
a further quantum speedup is possible (red solid line). 



Our proposed strategy is to construct a coherent encod- 
ings of a thermal state (GETS) \ipcETs) directly, rather 
than sampling from the thermal probability distribution: 



1^, 



CETS) 



) = Y,^e-P"(^)/Z\s). 



(1) 



where s = {0,1}^, f3 = l/ksT is the inverse temper- 
ature, H{s) is the eigen-energy of some classical spin 
Hamiltonian for the A/'-spin configuration s = siS2 ■ - ■ sn 
and Z is the partition function. This CETS can be 
transformed into the corresponding thermal state p = 
e~^^/Tr (e~^^) by including a set of N ancilla qubits, 
performing bit-by-bit GNOTs such that |s) (g) |0 • • • 0)^ ^ 



k) 



and tracing over the ancilla system. However, 



for some applications, such as the partition functions es- 
timation in [11 it is preferable to use the GETS directly. 
Below we present a method for preparing the GETS of 
a classical Hamiltonian from the initial state |0 • • • 0) by a 
sequence of locally-controlled rotations. Zalka's approach 
[To] , as applied to discrete cases [12] , allows for preparing 
the GETS by adding qubits one by one, and performing 
a rotation (controlled by all of the previous qubits) on 
each new qubit as 

ki • • • Sk) |0) ^ |5i • • • Sk) {cosOs |0) + sin^, |1)) , (2) 

where cos^ Os is the conditional probability of Sk-\-i = 0, 
given that the first k spins are in a particular configu- 
ration S1S2 ■ - ■ Sk. The problem here is that in general, 
this requires the knowledge (or efficient calculation) of 
O (2^) conditional probabilities. Thus, Zalka's method 
is efficient only when the probability distributions are 
efficiently integrable flT . Here we focus on the cases 
where the controlled rotations are local, i.e., they de- 
pend only on a few previous qubits. This in turn allows 
efficient computation of the respective rotation angles. 
More precsiely, we take into account the structure of the 
geometry of the physical systems and use the idea of the 
renormalization approach to obtain the rotation angles. 
Real-space renormalization — This method is also re- 
lated to the renormalization group method [M], which 
idea is to integrate out some degrees of freedom (coarse- 
graining) in the partition function Z, and describe the 
sub-system with a similar system with modiffed (renor- 
malized) couplings. As an example, consider a linear 
chain of three spins (Fig [2k). The partition function af- 
ter eliminating spin 3 (cf. Eq.([T3|)), 



A(/3)E 



,B(/3)siS2 



(3) 



is proportional to that of spin 1 and spin 2 interacting 
with an effective interaction —B (/3) //3. In contrast to 
this conventional renormalization treatment, where the 
degrees of freedom of the physical systems are progres- 
sively reduced, our method works in a reverse fashion: at 
each step, we increase the number of degrees of freedom, 
and then perform a controlled-rotation (Fig. [2b), which 
also changes the effective interaction of spins 1 and 2. 

Below, we define a sequential method for preparing a 
GETS for a generalized Ising Hamiltonian of N classical 
spins which has multiple spin-spin coupling constants: 



H. 



7 v^j^j 



E^. 



7 ^ ^ijk^i^j^k 
ijk 



(4) 



Our goal is to investigate how a GETS can be constructed 
by locally-controlled quantum rotations. Suppose we are 
given a GETS as defined in Eq. (fTl) of /c spins {ipk) foi" the 
Hamiltonian given by Eq. Q, and an additional qubit 
initialized in the state |0) that will become spin k -\- 1 of 
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FIG. 2: (color online) Real space renormalization approach 
for preparing the coherent encoding of a thermal state 
(GETS), (a) Spin 3 is eliminated by integration, inducing an 
effective interaction (blue bond) between spin 1 and spin 2. 
(b) A controlled rotation is performed on spin 3, inducing an 
effective interaction (red bond) between spin 1 and spin 2. (c) 
A quantum circuit demonstrating the sequential construction 
of the full thermal state. 



our GETS. Let us define the rotation angle, ^5, by 



-I3m, 



/Ws, 



(5) 



where Ws = e~^^^ + e^^^ = 2 cosh (/3ms), and rris = 
771 (51, . . . , 5/c) is a function of the spin variables of the 
first k spins. After performing a controlled rotation ^ 
on spin /c + 1, with angles given by Eq. ([5|, we obtain a 
GETS IV^/c+i) of a new {k + l)-spin Hamiltonian 

^fc+i = Hk^rris- Sfc+i. (6) 

To justify this statement, rewrite W in Eq. ^ as 



Ws 



-(3ms 



of^^s 



,^-p{H,-H,) 



(7) 



for some constant A^ and some fc-spin Ising spin Hamil- 
tonian Hk (with possible higher-order interactions). The 
state that we get from |^fe) |0) by the controlled rotation 
(d) is Y.sVPts)\s) with F{s) = F(si,S2,...,s„) given 

by 



F{s) 



^-(3Hk g-/3ms-Sfc+i 



w:. 



-l3(^Hk-\-ms-Sk+i) 
Zk^k 



(8) 



i.e. a GETS for the Hamiltonian (Ig]). Moreover, the 
new normalization constant is the same as the partition 
function Z^j^i for the system with Hamiltonian ^ and 
can be obtained simply by Z^+i = A^Zk- The term 
Hk in ^ is an Ising Hamiltonian of the form Q for 
the first k spins, but associated with a different set of 
renormalized couplings {Aj^ Bij^Cijk^ •'•}' Finally, the 



constant A^ can be shown to be the geometric mean of 
the LHS of ^ 



A/e = 2 ]^ [cosh {2Pms 



1I/2" 



(9) 



This is reminiscent of formulas which appear in classical 
algorithms such as belief propagation [8] for calculating 
some thermal properties of some spin systems. The con- 
trolled rotation is therefore the crucial element of our 
renormalization step. Using our method iteratively as 
shown in Fig. [2J3, we can generate the GETS of a partic- 
ular spin Hamiltonian. 

In general, Hk could contain up to /c-local interaction 
terms. If all the terms in Hk involve at most t spins, 
we call this a t-renormalizable operation. From now on, 
we will restrict ourselves to 2-renormalizable operations, 
which for example, includes local magnetic fields and 
two-spin interactions. Next, we will consider a general 
construction of a GETS for systems with finite range in- 
teraction, followed by some specific examples for further 
illustration. 

Finite-range interactions and belief propagation — As 
a general construction, we consider spin chains with 
finite-range interactions involving z neighboring spins. 
The computational complexity of this approach gener- 
ally scales exponentially in z. As an example, consider 
two groups of spins s and t, each can be considered as a 
2^ dimensional system. The Hamiltonian is of the form 



H = Hs -\- Ht-\- Hs 



(10) 



where Hs and Ht are the internal interaction terms for 
spins within group s and t, and Hgt contains the inter- 
actions between the groups. We start with preparing the 
state of the group s as M~^/^ ^^ ye^^^^T^ |s), where 
M is a normalization constant, and 7^ is a function of 
the spins s to eliminate renormalization effects induced 
by the spins in group t. The group t is initialized in the 
state |0 • • • 0). We choose the controlled rotation Q, 



\s) |0 . . . 0) ^ |5) ^ ^e-«^'+«-) • (7t/7«) \t) , (11) 

t 

with jt determined by the next group of spins to be in- 
cluded in the preparation procedure. If group t is the 
last group, then all jt are equal to 1. To ensure unitarity 
of this operation, we require 



^7,e-'5(^*+^-), 



(12) 



which is a recursion relation typically encountered in be- 
lief propagation [8 problems. For a group of z spins, and 
a given set of 7^, the sum involves 0(2^) terms, scaling 
exponentially in z. To perform the multi-qubit rotation, 
we can apply Zalka's algorithm [10], which requires the 
computation of 0(2^) rotation angles, and a polynomial 



number of subsequent quantum operations. To save com- 
putational resource for large z, it is more efficient to de- 
termine the angles for rotation "on the ffy" . This can be 
achieved by the quantum amplitude amplification algo- 
rithm [^ calculated with some ancilla qubits. 

We can apply this approach to a^n N x N square lattice 
of Ising spins with non-uniform couplings and arbitrary 
local magnetic fields. We make a group for each row of 
z = N spins. In the worst case scenario, the number of 
required operations in the above approach then scales [18] 
as 0(2^^), which becomes 0(2^) after combining with 
the amplitude amplification algorithm. This is still an ex- 
ponential algorithm, but nevertheless with an exponen- 
tial speed up over the direct application of Zalka's algo- 
rithm, whose complexity scales as 0(2^ ), as it requires 
the preparation of a probability distribution with 2^ 
amplitudes. However, for the uniform 2D Ising model 
without magnetic fields, an efficient t-renormalizable ap- 
proach might exist, as classical polynomial algorithms 
exist for this problem [16 . 

In the following, we illustrate our method by explicitly 
giving several examples of physical interest. 

Building blocks for frustrated magnets and spin ice 
— As the first example, we show how to generate a 
GETS of a triangle plaquette of three Ising spins by a 
2-renormalizable operation. Our goal is to prepare a 
GETS of three spins (see Fig. |2]3), for the Hamiltonian 
Hs = JsiS2-\-JsiSs-\-Js2Ss. Let us start with two qubits 
initialized as: M'^ Esi,s2={o,i} VlsiS2e-^'^'''^ \s1S2), 
where M is a normalization constant and jsiS2 > is 
some positive function of si and S2 to be determined 
later. Let us add a third qubit in the state |0) to the sys- 
tem, and act with the controlled rotation ([5| depending 
on the values of the first two qubits. When we choose 
rus = J{si + S2) for some constant J, we can use the 
well-known result in renormalizing the ID Ising chain 
[H], and write 



Ws 



-/3J(si+s2) ^ g^J(si+s2) ^ Ae/^^^ 



where the coefficients A and B are 



A 
B 



2v'cosh(2/3J), 
(1/2/3) In cosh (2/3 J). 



(13) 



(14) 



Observe now that if we chose jsis 



^^(3Bs^S2 ^l^Q^ 



preparing the first two qubits, applying the controlled 
rotation of the third qubit eliminates this factor. Gonse- 
quently, this operation produces the GETS for the 3-spin 
Ising cycle Hs. 

Furthermore, we can build on top of this triangular 
lattice to prepare a GETS of the basic tetrahedral unit 
of spin ice [17 (see Fig. [3]^ a)). Let's add a fourth qubit 
to the system, and perform a rotation controlled by the 
first three qubits: 

\s1s2ss) |0) ^ \s1s2ss) {cosOs |0) + sin^, |1)) , (15) 




FIG. 3: (color online) Examples of thermal state generation 
for the basic units of Ising spin cycles and spin ice (a) A 
rotation on spin 4, controlled by spin 1,2, and 3, results a 
thermal state of the tetrahedral spin ice. (b). A next- nearest- 
neighbor spin chain is created by consecutive application of 
two-qubit controlled rotations. 



where cos Os is given in ^ with 
In this case, we have, 



J{si 



S2 + S3). 



W.. 



^^K{siS2+S2S3-\-SiS3) 



(16) 



^1/4 



A = 2(cosh^(J)cosh(3J))' 

K = (1/4) In (cosh (3J)/ cosh (J)). 



This suggests that to compensate for the renormalization 
effect of adding the fourth spin, we can replace J with 



J - K in ( 13 ) and (h4l when preparing the GETS of the 



base triangle of the spin-ice tetrahedron that is indicated 
as red bonds in Fig. IsFa). 

Ising spin chains with next-nearest-neighbor interac- 
tions and local magnetic fields — As another application 
of our renormalization method, we consider the applica- 
tion to spin chain with next-nearest-neighbor interaction, 
which generalizes the work done in [9 . Given a chain 
of n spins with nearest-neighbor interactions J^, next- 
nearest-neighbor interactions L^, and local fields /i^, the 
Hamiltonian is, 



n-l 



n-2 



H = Y^ JiSiSi^i + ^ LiSiSi^2 + y^ hiSi. 









(17) 



For i > 3, we define a 2-renormalizable controlled- 
rotation on the i-th qubit: 

\si-2Si-1) |0) ^ \s^-2Si-i) (cos^, |0) + sin^, |1)) , (18) 

with rris = Li-2Si-2 + Ji-iSi-i + hiSi in ([s]) determining 
the rotation angles. Using 

gf^ = 2 cosh {hi ± J,_i ± L,_2) , (19) 

and following the procedure in ([7|), we can rewrite 



This equality is satisfied for the following choices 

A. = {grgt-g-^grf, (21) 

Ci = {l/4P)Hg+-g++/g--g-+), 
A = (1/4/3) In (5-+5++/5--5+-). 

After the controlled rotation on the i-th qubit, the spin- 
spin coupling between the spins z — 2, z — 1 is thus renor- 
malized as Ji-2 -^ Ji-2 + Bi. Also, the local magnetic 
fields for these spins are renormalized as hi-i -^ hi-i-\-Ci 
and hi-2 -^ hi-2 + A- Note that the next-nearest- 
neighbor coupling Li_2 term remains unchanged. 



Our goal is a GETS of the Hamiltonian (17). We can 



eliminate the unwanted renormalization effects with the 
following procedure: (i) for i = l..n, initialize h[i] = hi, 
J[i] = Ji and L[i] = L^; (ii) compute Bi, Ci and Di 
from (21) using the values stored in (h[i], J[i], L[i]); put 
h[i - 1] = h[i - 1] - Ci, h[i - 2] = h[i - 2] - A and 
J[i — 2] = J[± — 2] — Bi] let i = i — 1. (iii) Repeat step 
(ii) until i = 3. The desired GETS can then be obtained 
by a sequence of 2-locally controlled rotations using the 
stored values in (h[i], J[i], L[i]) instead of the original 
hi , Ji and Li . 

Conclusion — To summarize, we have developed an 
algorithm which identifies a class of classical spin prob- 
lems that can be simulated efficiently with a quantum 
computer. In this class of problems, our method scales 
efficiently compared with MGMG methods, as it is inde- 
pendent of the gap of the Markov chain and temperature. 
On the other hand, we believe that the tools developed 
here could be useful for classifying the complexity classes 
of certain spin models. An avenue for further research 
is the complexity classification of spin systems by their 
t-renormalizability, which may suggest a deeper under- 
standing of the connection between complexity theory 
and quantum simulation. 
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